#include <iostream>
#include "spline.h"
#include "matplotlibcpp.h"
#include <math.h>

namespace plt = matplotlibcpp;

double func(double x)
{
    return 1 / (1+25*x*x);
}

double dfunc(double x)
{
    return (-50 * x) / ((1+25*x*x)*(1+25*x*x));
}

double ddfunc(double x)
{
    return 50*(75*x*x-1) / (pow((1+25*x*x),3));
}

int main()
{

    double a = -1.0, b = 1.0;
    std::vector<int> N = {6,11,21,41,81} ;
    //not_a_knot
    std::vector<spline> S1;
    S1.resize(N.size());
    for (int k = 0; k < N.size(); k++)
    {
        std::vector<double> x(N[k]);
        std::vector<double> f(N[k]);
        double gap = (b - a)/(N[k] - 1);
        for(int i = 0; i < N[k]; i++)
        {
            x[i] = a + i * gap;
            f[i] = func(x[i]);
        }
        spline temp(x,f);
        S1[k] = temp;
        S1[k]. not_a_knot();

    }
    
    //以下为画图部分
    //func
    int n = 1000;    //描点画图的点数
    double figgap = (b - a) / n;
    std::vector<double> u(n+1), y(n+1);

    plt::figure();          //创建一张图

    for (int k = 0; k < N.size(); k++)
    {
        for (int j = 0; j <= n; j++)
        {
            u[j] = a + j * figgap;
            y[j] = S1[k].fit(u[j]);
        }
        plt::plot(u, y,{{"label", "not_a_knot N="+to_string(N[k])}, {"ls", "--"}});                //
    }
    
    for (int j = 0; j <= n; j++)
    {
        u[j] = a + j * figgap;
        y[j] = func(u[j]);
    }
    plt::plot(u, y,{{"label", "$y = \\frac{1}{1+25x^2} $"},{"color","k"}});
        
    plt::xlim(a, b);
    plt::title("not_a_knot spline in A"); // set a title
    plt::legend();                // enable the legend

    plt::save("../fig/A-not_a_knot-func.png");
    std::cout << "output : A-not_a_knot-func.png"  << endl;

    //error
    plt::figure();          //创建一张图

    for (int k = 0; k < N.size(); k++)
    {
        for (int j = 0; j <= n; j++)
        {
            u[j] = a + j * figgap;
            y[j] = S1[k].fit(u[j])-func(u[j]);
        }
        plt::plot(u, y,{{"label", "not_a_knot N="+to_string(N[k])}, {"ls", "--"}});                //
    }
        
    plt::xlim(a, b);
    plt::title("not_a_knot spline error in A"); // set a title
    plt::legend();                // enable the legend

    plt::save("../fig/A-not_a_knot-error.png");
    std::cout << "output : A-not_a_knot-error.png"  << endl;

    //convergence rate
    plt::figure();          //创建一张图

    std::vector<double> errs, num_of_inter;

    for (int k = 0; k < N.size(); k++)
    {
        num_of_inter.push_back(log(N[k]-1));
        double err = 0;
        double gap = (b - a)/(N[k] - 1);
        for (int j = 0; j < N[k]-1; j++)
        {
            double temp = fabs(S1[k].fit(a+j*gap+gap/2) - func(a+j*gap+gap/2));
            if (temp > err) err = temp;
        }                //
        errs.push_back(log(err));
    }
        
    plt::plot(num_of_inter,errs,{{"label", "ln(err) - ln(num_of_inter)"}});

    plt::title("not_a_knot spline convergence rate in A"); // set a title
    plt::legend();
    plt::save("../fig/A-not_a_knot-rate.png");
    std::cout << "output : A-not_a_knot-rate.png"  << endl;

}